from utopia.monte_carlo.montecarlo_analysis_paper import run_mc_analysis_paper, set_emission, sim_to_dataframe, reshape_mc_results
from utopia.monte_carlo.monte_carlo_helpers import *
from utopia.utopia import utopiaModel
import numpy as np
import scipy.stats as st
import pandas as pd

c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\monaco\dvars_sensitivity.py:18: UserWarning: Consider installing numba for better performance
  warn("Consider installing numba for better performance", UserWarning)

data_data = utopiaModel.load_json_file("data/default_data.json")
config_data = utopiaModel.load_json_file("data/default_config.json")

import numpy as np
# To reproduce Rakesh results, use these following discrete values to then generate distributions (lognormal and uniform) from which to sample:

# Plastic density (mpdensity): [900, 1100, 1500] kg/m³ 
# Fragmentation style (FI): [0, 0.3, 0.5, 0.7, 1.0] - ranging from more erosive to more explosive fragmentation 
# Discorporation half time (thalf_deg): [316,1000,10000, 31600, 100000] days 
# Fragmentation timescale (thalf_frag): [10, 31.6, 100, 316] days
# your discrete support values

vals_density = np.array([900, 1100, 1500])
vals_FI = np.array([0, 0.3, 0.5, 0.7, 1.0])
vals_deg = np.array([316, 1000, 10000, 31600, 100000])
vals_frag = np.array([10, 31.6, 100, 316])

def lognorm_params_from_values(values, sigma):
    """
    Return (s, scale) for scipy.stats.lognorm given target values.
    sigma = shape parameter (spread in log space).
    We center the lognormal at the geometric mean of the values.
    """
    mu = np.mean(np.log(values))  # geometric mean in log-space
    scale = np.exp(mu - (sigma**2) / 2)
    return sigma, scale

# choose sigma to control spread (tune manually!)
s_deg, scale_deg = lognorm_params_from_values(vals_deg, sigma=0.5)
s_frag, scale_frag = lognorm_params_from_values(vals_frag, sigma=0.8)
s_density, scale_density = lognorm_params_from_values(vals_density, sigma=0.1)
s_FI, scale_FI = lognorm_params_from_values(vals_FI[vals_FI > 0], sigma=0.5)  # avoid zero for lognorm

param_distributions = {
    "MPdensity_kg_m3": ("lognorm", {"s": s_density, "scale": scale_density}),
    "FI": ("uniform", {"loc": 0.0, "scale": 1.0}),
    "t_half_deg_free": ("lognorm", {"s": s_deg, "scale": scale_deg}),
    "t_frag_gen_FreeSurfaceWater": ("lognorm", {"s": s_frag, "scale": scale_frag})
}

# Visualise the distributions by plotting

for psrsm in param_distributions:
    plot_distribution_from_dict(psrsm, param_distributions)

# Recent studies have challenged the density of biofilm used here (1388 kg m-3) as being too high. see e.g. # https://doi.org/10.1186/s43591-024-00098-2 
# Consider testing chnages in biofouling density

# Set emission scenarios to run that are compatible for comparison with observations data from Kooi et al.

##############################################
# DEFINE ARTICLE -> EMISSION SCENARIOS
##############################################

overlay_emission_scenarios = {
    "Sea Sediment | Bergmann":["Ocean_Surface_Water","Coast_Surface_Water","Coast_Column_Water"],
    "Sea Surface all | Cai":["Surface_Freshwater","Coast_Surface_Water","Impacted_Soil_Surface"],
    "Sea Surface Fragments | Enders":["Ocean_Surface_Water","Coast_Surface_Water"],
    "Sea Surface Fibers | Enders":["Ocean_Surface_Water","Coast_Surface_Water"],
    "Fragments in water | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Fibers in water | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Fragments in sediment | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Fibers in sediment | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Sea Surface all | Erni-Cassola":["Coast_Surface_Water","Impacted_Soil_Surface"],
    "Fresh water sediments | Imhof":["Surface_Freshwater","Impacted_Soil_Surface"],
    "Sea Surface all | Isobe L1":["Surface_Freshwater","Coast_Surface_Water"],
    "Sea Surface all | Isobe L2":["Surface_Freshwater","Coast_Surface_Water"],
    "Sea Surface all | Isobe L3":["Surface_Freshwater","Coast_Surface_Water"],
    "Sea Surface all | Isobe L4":["Surface_Freshwater","Coast_Surface_Water"],
    "FP all | Scheurer":["Impacted_Soil_Surface"],
    "Sea Surface Fragments | Song":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Sea Surface Fibers | Song":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Sea Surface all | Zhang":["Coast_Surface_Water","Surface_Freshwater"]}
    
    
    
    
# DEFINE ARTICLE->COMPARTMENTS MAPPING
#############################################
overlay_compartments_mapping = {
    "Sea Sediment | Bergmann": ["Sediment_Coast", "Sediment_Ocean"],
    "Sea Surface all | Cai": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface Fragments | Enders": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface Fibers | Enders": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Fragments in water | Eo": ["Bulk_Freshwater", "Surface_Freshwater"],
    "Fibers in water | Eo": ["Bulk_Freshwater", "Surface_Freshwater"],
    "Fragments in sediment | Eo": ["Sediment_Ocean", "Sediment_Coast", "Beaches_Deep_Soil","Beaches_Soil_Surface", "Sediment_Freshwater"],
    "Fibers in sediment | Eo": ["Sediment_Ocean", "Sediment_Coast", "Beaches_Deep_Soil",
    "Beaches_Soil_Surface", "Sediment_Freshwater"],
    "Sea Surface all | Erni-Cassola": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Fresh water sediments | Imhof": ["Bulk_Freshwater", "Surface_Freshwater"],
    "Sea Surface all | Isobe L1": ["Coast_Surface_Water", "Ocean_Surface_Water"], 
    "Sea Surface all | Isobe L2": ["Coast_Surface_Water", "Ocean_Surface_Water"], 
    "Sea Surface all | Isobe L3": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface all | Isobe L4": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "FP all | Scheurer": ["Coast_Surface_Water", "Ocean_Surface_Water"], 
    "Sea Surface Fragments | Song": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface Fibers | Song": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface all | Zhang": ["Coast_Surface_Water", "Ocean_Surface_Water"]
}

# NOTE! Maybe the mapping for the FP all | Scheurer study has to be revisited. Samples were collected in floodplain soils, so the impacted_soil_surface compartment might be more appropriate. Since the model does not simulate floodplain soils as a separate compartment.?

# Set emission scenarios according to the mapping defined above to the run the model for each parameter combination using the Monte Carlo framework and perform the comparison with the observed data according to mapping.

import seaborn as sns
import matplotlib.pyplot as plt
from utopia.monte_carlo.montecarlo_analysis_paper import run_and_plot_top_results

MC_simulations_results_all={}
for art in overlay_emission_scenarios.keys():
    print(f"Running Monte Carlo analysis for article: {art}")
    MC_simulations_results=[]
    for comp in overlay_emission_scenarios[art]:
        # Set emissions to that compartment and include it in base_data file
        data_data["emiss_dict_g_s"]= set_emission(data_data["emiss_dict_g_s"].copy(), comp, "e", 100)
        
        
        #Run Monte Carlo analysis for each emission scenario and targer compartment
        for comp_T in overlay_compartments_mapping[art]:
            emission_comp=comp
            target_comp=comp_T
            sim= run_mc_analysis_paper(
            base_config=config_data,
            base_data=data_data,
            param_distributions=param_distributions,
            n_cases=20,art=art,emission_comp=emission_comp,target_comp=target_comp
        )
        
            # Convert sim results to dataframe
            df_MC_results = sim_to_dataframe(sim)
            
            
            MC_simulations_results.append(df_MC_results)
    MC_simulations_results = pd.concat(MC_simulations_results, ignore_index=True)
    input_params=["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]
    tidy_mc_results=reshape_mc_results(MC_simulations_results, input_params, art)
    results_sorted_art = tidy_mc_results.sort_values("RMSE", ascending=True)
    # Keep e.g. top 20 best-fitting combinations
    top_results_art = results_sorted_art.head(20)
    run_and_plot_top_results(top_results_art[0:2],data_data, config_data)
    

    # Correlation heatmap for the top results
    corr = top_results_art[0:10][["RMSE", "R2", "Slope", "MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]].corr()

    plt.figure(figsize=(8,6))
    sns.heatmap(corr, annot=True, cmap="coolwarm", center=0)
    plt.title("Correlation matrix (Top 10 best fits)")
    plt.show()
    
    
    MC_simulations_results_all[art]=top_results_art

Running Monte Carlo analysis for article: Sea Sediment | Bergmann
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 4997.68 cases/s]
Running cases:  30%|███       | 6/20 [00:13<00:31,  2.22s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:19<00:23,  2.18s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:44<00:00,  2.21s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 1392.18 cases/s]

Simulation complete! Runtime: 0:00:44.245699
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<?, ? cases/s]
Running cases:  30%|███       | 6/20 [00:13<00:31,  2.22s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:19<00:24,  2.25s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:45<00:00,  2.26s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 12685.03 cases/s]

Simulation complete! Runtime: 0:00:45.289060
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 5318.67 cases/s]
Running cases:  30%|███       | 6/20 [00:13<00:31,  2.28s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:20<00:25,  2.28s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:47<00:00,  2.38s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<?, ? cases/s]

Simulation complete! Runtime: 0:00:47.558976
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 13856.31 cases/s]
Running cases:  30%|███       | 6/20 [00:17<00:41,  2.98s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:25<00:31,  2.86s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:49<00:00,  2.49s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<?, ? cases/s]

Simulation complete! Runtime: 0:00:49.858126
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 5517.01 cases/s]
Running cases:  30%|███       | 6/20 [00:13<00:31,  2.22s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:19<00:24,  2.19s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:45<00:00,  2.26s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 1077.66 cases/s]

combined_df = pd.concat(MC_simulations_results_all.values(), ignore_index=True)

combined_df

combined_results_shorted=combined_df.sort_values("RMSE", ascending=True)
combined_results_shorted

from utopia.monte_carlo.montecarlo_analysis_paper import run_and_plot_top_results
run_and_plot_top_results(combined_results_shorted[0:2],data_data, config_data)

import seaborn as sns
import matplotlib.pyplot as plt

# Correlation heatmap for the top results
corr = combined_results_shorted[0:10][["RMSE", "MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]].corr()

plt.figure(figsize=(8,6))
sns.heatmap(corr, annot=True, cmap="coolwarm", center=0)
plt.title("Correlation matrix (Top 10 best fits)")
plt.show()

# Pair-wise relationships for the top results
sns.pairplot(
    combined_results_shorted[0:10],
    vars=["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater","RMSE"],
    hue="Target_Compartment",
    corner=True
)
plt.suptitle("Parameter relationships for best-fitting runs", y=1.02)
plt.show()

# Read data from CSV if needed
combined_results = pd.read_csv("loop_all.csv")
df=combined_results.copy()

combined_results_shorted=df.sort_values("RMSE", ascending=True)

combined_results_shorted[1:20]

import seaborn as sns
import matplotlib.pyplot as plt

# Correlation heatmap for the top results
corr = combined_results_shorted[1:20][["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]].corr()

plt.figure(figsize=(8,6))
sns.heatmap(corr, annot=True, cmap="coolwarm", center=0)
plt.title("Correlation matrix (Top 20 best fits)")
plt.show()

# Pearson & Spearman correlation between Density and FI for all results

# Pearson tests linear correlation.

# Spearman tests monotonic correlation (rank-based, robust to nonlinearities).

# If p < 0.05, the correlation is statistically significant.

from scipy.stats import pearsonr, spearmanr

x = df["MPdensity_kg_m3"]
y = df["FI"]

pearson_corr, pearson_p = pearsonr(x, y)
spearman_corr, spearman_p = spearmanr(x, y)

print(f"Pearson r = {pearson_corr:.3f}, p = {pearson_p:.3e}")
print(f"Spearman ρ = {spearman_corr:.3f}, p = {spearman_p:.3e}")

import seaborn as sns
import matplotlib.pyplot as plt

sns.lmplot(
    data=df[1:20], 
    x="MPdensity_kg_m3", 
    y="FI", 
    #hue="Target_Compartment", 
    height=6, aspect=1.2,
    scatter_kws={"alpha": 0.7}
)
plt.title("Density vs FI (with regression fits)")
plt.show()

Pearson r = 0.140, p = 7.955e-03
Spearman ρ = 0.016, p = 7.579e-01

from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.model_selection import train_test_split

# Suppose your dataframe is df
# Identify features and target
target = 'RMSE'
numeric_features = ['MPdensity_kg_m3', 'FI', 't_half_deg_free', 't_frag_gen_FreeSurfaceWater']
#categorical_features = ['Observed_dataset', 'Emission_Compartment', 'Target_Compartment']

#X = df[numeric_features + categorical_features]
X = df[numeric_features]
y = df[target]

# One-hot encode categorical columns
preprocessor = ColumnTransformer(
    transformers=[
        ('num', 'passthrough', numeric_features),
        #('cat', OneHotEncoder(handle_unknown='ignore'), categorical_features)
    ]
)

from sklearn.ensemble import RandomForestRegressor
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt
import pandas as pd

# Build pipeline
model = Pipeline(steps=[
    ('preprocess', preprocessor),
    ('rf', RandomForestRegressor(n_estimators=500, random_state=42))
])

# Train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# Fit model
model.fit(X_train, y_train)

# Extract feature importances
rf = model.named_steps['rf']
feature_names = model.named_steps['preprocess'].get_feature_names_out()
importances = pd.Series(rf.feature_importances_, index=feature_names).sort_values(ascending=False)

# Plot importance
plt.figure(figsize=(10, 5))
importances[:15].plot(kind='barh')
plt.title('Feature Importance (Random Forest)')
plt.xlabel('Importance')
plt.gca().invert_yaxis()
plt.show()

importances.head(10)

num__FI                             0.373965
num__t_half_deg_free                0.247732
num__MPdensity_kg_m3                0.190985
num__t_frag_gen_FreeSurfaceWater    0.187317
dtype: float64

import seaborn as sns
corr = df[numeric_features + [target]].corr()
plt.figure(figsize=(8, 6))
sns.heatmap(corr, annot=True, cmap='coolwarm', fmt=".2f")
plt.title("Correlation matrix (numeric variables)")
plt.show()

corr[target].sort_values(ascending=False)

RMSE                           1.000000
MPdensity_kg_m3                0.237672
t_frag_gen_FreeSurfaceWater    0.155295
t_half_deg_free                0.091898
FI                             0.075080
Name: RMSE, dtype: float64

from utopia.monte_carlo.montecarlo_analysis_paper import run_and_plot_top_results
run_and_plot_top_results(top_results_all[0:5],data_data, config_data)

import seaborn as sns
import matplotlib.pyplot as plt

# Correlation heatmap for the top results
corr = top_results[["RMSE", "R2", "Slope", "MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]].corr()

plt.figure(figsize=(8,6))
sns.heatmap(corr, annot=True, cmap="coolwarm", center=0)
plt.title("Correlation matrix (Top 10 best fits)")
plt.show()

# Pair-wise relationships for the top results
sns.pairplot(
    top_results,
    vars=["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater","RMSE"],
    hue="Target_Compartment",
    corner=True
)
plt.suptitle("Parameter relationships for best-fitting runs", y=1.02)
plt.show()

fig, axes = plt.subplots(2, 2, figsize=(10, 8))

params = ["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]

for ax, param in zip(axes.ravel(), params):
    sns.scatterplot(data=top_results, x=param, y="RMSE", hue="Target_Compartment", ax=ax)
    ax.set_title(f"RMSE vs {param}")

plt.tight_layout()
plt.show()

# Pearson & Spearman correlation between Density and FI for all results

# Pearson tests linear correlation.

# Spearman tests monotonic correlation (rank-based, robust to nonlinearities).

# If p < 0.05, the correlation is statistically significant.

from scipy.stats import pearsonr, spearmanr

x = top_results["MPdensity_kg_m3"]
y = top_results["FI"]

pearson_corr, pearson_p = pearsonr(x, y)
spearman_corr, spearman_p = spearmanr(x, y)

print(f"Pearson r = {pearson_corr:.3f}, p = {pearson_p:.3e}")
print(f"Spearman ρ = {spearman_corr:.3f}, p = {spearman_p:.3e}")

import seaborn as sns
import matplotlib.pyplot as plt

sns.lmplot(
    data=top_results, 
    x="MPdensity_kg_m3", 
    y="FI", 
    #hue="Target_Compartment", 
    height=6, aspect=1.2,
    scatter_kws={"alpha": 0.7}
)
plt.title("Density vs FI (with regression fits)")
plt.show()

Pearson r = 0.114, p = 6.323e-01
Spearman ρ = 0.130, p = 5.844e-01

# Pearson & Spearman correlation between Density and FI for best fits (e.g. lowest 10% RMSE)

threshold = top_results["RMSE"].quantile(0.1)  # lowest 10% of RMSE
best_fits = top_results[top_results["RMSE"] <= threshold]

pearson_corr, _ = pearsonr(best_fits["MPdensity_kg_m3"], best_fits["FI"])
print(f"Best fits Pearson r = {pearson_corr:.3f}")


sns.lmplot(
    data=best_fits, 
    x="MPdensity_kg_m3", 
    y="FI", 
    hue="Target_Compartment", 
    height=6, aspect=1.2,
    scatter_kws={"alpha": 0.7}
)
plt.title("Density vs FI (with regression fits)")
plt.show()

Best fits Pearson r = 0.114

from scipy.stats import pearsonr, spearmanr
import pandas as pd

# Keep only best 10% RMSE across all compartments
threshold = top_results["RMSE"].quantile(0.1)
best_fits = top_results[top_results["RMSE"] <= threshold]

results = []

for comp in best_fits["Target_Compartment"].unique():
    subset = best_fits[best_fits["Target_Compartment"] == comp]
    
    if len(subset) > 2:  # need at least 3 points
        pearson_corr, pearson_p = pearsonr(subset["MPdensity_kg_m3"], subset["FI"])
        spearman_corr, spearman_p = spearmanr(subset["MPdensity_kg_m3"], subset["FI"])
    else:
        pearson_corr, pearson_p, spearman_corr, spearman_p = [None]*4
    
    results.append({
        "Target_Compartment": comp,
        "Pearson_r": pearson_corr,
        "Pearson_p": pearson_p,
        "Spearman_rho": spearman_corr,
        "Spearman_p": spearman_p,
        "n_points": len(subset)
    })

correlation_best = pd.DataFrame(results)
print(correlation_best)

    Target_Compartment Pearson_r Pearson_p Spearman_rho Spearman_p  n_points
0  Coast_Surface_Water      None      None         None       None         1
1      Bulk_Freshwater      None      None         None       None         1

import seaborn as sns
import matplotlib.pyplot as plt

# --- pick best fits per compartment (lowest 10% RMSE) ---
best_fits_list = []

for comp in top_results["Target_Compartment"].unique():
    subset = top_results[top_results["Target_Compartment"] == comp]
    if len(subset) > 2:
        threshold = subset["RMSE"].quantile(0.1)   # lowest 10%
        best_subset = subset[subset["RMSE"] <= threshold]
        best_fits_list.append(best_subset)

best_fits = pd.concat(best_fits_list, ignore_index=True)

# --- plot each compartment separately ---
g = sns.lmplot(
    data=best_fits,
    x="MPdensity_kg_m3",
    y="FI",
    col="Target_Compartment",   # separate plot per compartment
    col_wrap=2,                 # wrap panels (2 per row)
    height=4,
    scatter_kws={"alpha": 0.7},
    line_kws={"color": "black"}
)

g.set_titles("{col_name}")
plt.suptitle("Density vs FI (Best RMSE fits, per compartment)", y=1.05)
plt.show()

sim.outvars

{'Spearman_r': <monaco.mc_var.OutVar at 0x16cc2281650>,
 'Spearman_p': <monaco.mc_var.OutVar at 0x16cbfaadc90>,
 'Pass_Spearman': <monaco.mc_var.OutVar at 0x16cbd9294d0>,
 'RMSE': <monaco.mc_var.OutVar at 0x16cbfb8df50>,
 'R_squared': <monaco.mc_var.OutVar at 0x16cbfb8c990>,
 'Model_slope': <monaco.mc_var.OutVar at 0x16cc1deb590>,
 'Model_intercept': <monaco.mc_var.OutVar at 0x16cc1deba90>,
 'Target_compartment': <monaco.mc_var.OutVar at 0x16cc1dea490>,
 'Emission_compartment': <monaco.mc_var.OutVar at 0x16cc1de8550>}

sim.calcSensitivities("C_g_m3_Ocean_Surface_Water")
fig, ax = sim.outvars["C_g_m3_Ocean_Surface_Water"].plotSensitivities()

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[24], line 1
----> 1 sim.calcSensitivities("C_g_m3_Ocean_Surface_Water")
      2 fig, ax = sim.outvars["C_g_m3_Ocean_Surface_Water"].plotSensitivities()

File c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\monaco\mc_sim.py:916, in Sim.calcSensitivities(self, outvarnames, tol, verbose)
    913 outvarnames = get_list(outvarnames)
    915 for outvarname in outvarnames:
--> 916     if not self.outvars[outvarname].isscalar:
    917         vwarn(self.verbose, f"Output variable '{outvarname}' is not scalar," +
    918                              'skipping sensitivity calculations.')
    919     else:

KeyError: 'C_g_m3_Ocean_Surface_Water'

sim.calcSensitivities("persistence_number")
fig, ax = sim.outvars["persistence_number"].plotSensitivities()

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[25], line 1
----> 1 sim.calcSensitivities("persistence_number")
      2 fig, ax = sim.outvars["persistence_number"].plotSensitivities()

File c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\monaco\mc_sim.py:916, in Sim.calcSensitivities(self, outvarnames, tol, verbose)
    913 outvarnames = get_list(outvarnames)
    915 for outvarname in outvarnames:
--> 916     if not self.outvars[outvarname].isscalar:
    917         vwarn(self.verbose, f"Output variable '{outvarname}' is not scalar," +
    918                              'skipping sensitivity calculations.')
    919     else:

KeyError: 'persistence_number'

import matplotlib.pyplot as plt

for var in sim.outvars:
    print(f"Variable: {var}")
    sim.calcSensitivities(var)
    fig, ax = sim.outvars[var].plotSensitivities()

    plt.show()

Variable: Spearman_r
Calculating sensitivity indices for 'Spearman_r'...
Calculating optimal hyperparameters Φ for 'Spearman_r' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

Variable: Spearman_p
Calculating sensitivity indices for 'Spearman_p'...
Calculating optimal hyperparameters Φ for 'Spearman_p' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

Variable: Pass_Spearman
Calculating sensitivity indices for 'Pass_Spearman'...
Calculating optimal hyperparameters Φ for 'Pass_Spearman' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

Variable: RMSE
Calculating sensitivity indices for 'RMSE'...
Calculating optimal hyperparameters Φ for 'RMSE' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

Variable: R_squared
Calculating sensitivity indices for 'R_squared'...
Calculating optimal hyperparameters Φ for 'R_squared' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

Variable: Model_slope
Calculating sensitivity indices for 'Model_slope'...
Calculating optimal hyperparameters Φ for 'Model_slope' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

Variable: Model_intercept
Calculating sensitivity indices for 'Model_intercept'...
Calculating optimal hyperparameters Φ for 'Model_intercept' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

Variable: Target_compartment
Calculating sensitivity indices for 'Target_compartment'...
Calculating optimal hyperparameters Φ for 'Target_compartment' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\monaco\dvars_sensitivity.py:277: RuntimeWarning: divide by zero encountered in log
  L = np.log(Rdet)/m + m*np.log(L_inner.T @ Rinv @ L_inner)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\scipy\optimize\_numdiff.py:686: RuntimeWarning: invalid value encountered in subtract
  df = [f_eval - f0 for f_eval in f_evals]
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\monaco\dvars_sensitivity.py:93: RuntimeWarning: invalid value encountered in divide
  ratios = sensitivities/sum(sensitivities)

Variable: Emission_compartment
Calculating sensitivity indices for 'Emission_compartment'...
Calculating optimal hyperparameters Φ for 'Emission_compartment' covariances...
Done calculating optimal hyperparameters.
Done calculating sensitivity indices.

sim.outvars.keys()

dict_keys(['Coast_Surface_Water_Spearman_r', 'Coast_Surface_Water_Spearman_p', 'Coast_Surface_Water_Pass_Spearman', 'Coast_Surface_Water_RMSE', 'Coast_Surface_Water_R_squared', 'Coast_Surface_Water_Model_slope', 'Coast_Surface_Water_Model_intercept', 'Ocean_Surface_Water_Spearman_r', 'Ocean_Surface_Water_Spearman_p', 'Ocean_Surface_Water_Pass_Spearman', 'Ocean_Surface_Water_RMSE', 'Ocean_Surface_Water_R_squared', 'Ocean_Surface_Water_Model_slope', 'Ocean_Surface_Water_Model_intercept', 'residence_time_mass', 'residence_time_number', 'persistence_mass', 'persistence_number'])

mc.multi_plot([sim.invars["MPdensity_kg_m3"],sim.outvars["persistence_mass"]],title="persistence_mass vs density",cov_plot=True)

(<Figure size 640x480 with 3 Axes>,
 (<Axes: xlabel='MPdensity_kg_m3'>,
  <Axes: ylabel='persistence_mass'>,
  <Axes: >))

	Case	Observed_dataset	Emission_Compartment	Target_Compartment	RMSE	R2	Slope	Spearman_r	Pass_Spearman	MPdensity_kg_m3	FI	t_half_deg_free	t_frag_gen_FreeSurfaceWater
0	108	Sea Sediment \| Bergmann	Coast_Column_Water	Sediment_Ocean	0.306841	-0.542379	-2.655127	-1.0	True	1334.859380	0.269348	2665.785177	31.280483
1	68	Sea Sediment \| Bergmann	Coast_Surface_Water	Sediment_Ocean	0.314419	-0.619506	-2.627525	-1.0	True	1334.859380	0.269348	2665.785177	31.280483
2	100	Sea Sediment \| Bergmann	Coast_Column_Water	Sediment_Ocean	0.319249	-0.669642	-2.453298	-1.0	True	1223.294660	0.054891	6510.335737	18.298327
3	112	Sea Sediment \| Bergmann	Coast_Column_Water	Sediment_Ocean	0.336169	-0.851313	-2.519515	-1.0	True	1170.689441	0.163532	10975.160488	41.971349
4	72	Sea Sediment \| Bergmann	Coast_Surface_Water	Sediment_Ocean	0.342700	-0.923945	-2.492947	-1.0	True	1170.689441	0.163532	10975.160488	41.971349
...	...	...	...	...	...	...	...	...	...	...	...	...	...
355	61	Sea Surface all \| Zhang	Surface_Freshwater	Ocean_Surface_Water	3.210394	-19.030002	-2.794831	-1.0	True	1027.531967	0.740818	4154.920691	47.156730
356	47	Sea Surface all \| Zhang	Surface_Freshwater	Coast_Surface_Water	3.554319	-23.551460	-2.733908	-1.0	True	1302.585757	0.539578	9481.499850	12.343862
357	23	Sea Surface all \| Zhang	Coast_Surface_Water	Ocean_Surface_Water	4.112485	-31.867982	-3.125205	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
358	63	Sea Surface all \| Zhang	Surface_Freshwater	Ocean_Surface_Water	4.395308	-36.544219	-3.319479	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
359	37	Sea Surface all \| Zhang	Coast_Surface_Water	Ocean_Surface_Water	4.651851	-41.054845	-3.557464	-1.0	True	1052.671922	0.811533	3073.415520	14.554329

	Case	Observed_dataset	Emission_Compartment	Target_Compartment	RMSE	R2	Slope	Spearman_r	Pass_Spearman	MPdensity_kg_m3	FI	t_half_deg_free	t_frag_gen_FreeSurfaceWater
120	149	Fragments in sediment \| Eo	Surface_Freshwater	Beaches_Deep_Soil	0.187272	0.865844	-1.988877	-1.0	True	1011.054694	0.966201	8317.906191	110.332414
220	6	Sea Surface all \| Isobe L2	Surface_Freshwater	Coast_Surface_Water	0.192740	0.540131	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
240	6	Sea Surface all \| Isobe L3	Surface_Freshwater	Coast_Surface_Water	0.196303	0.512105	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
121	137	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.210465	0.830557	-1.757989	-1.0	True	1052.671922	0.811533	3073.415520	14.554329
122	125	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.219214	0.816177	-1.605927	-1.0	True	1067.147583	0.814583	7649.971822	27.484251
...	...	...	...	...	...	...	...	...	...	...	...	...	...
295	34	FP all \| Scheurer	Impacted_Soil_Surface	Ocean_Surface_Water	6.005672	-108.220248	-4.254892	-1.0	True	1038.834733	0.450996	5968.470983	99.035129
296	32	FP all \| Scheurer	Impacted_Soil_Surface	Ocean_Surface_Water	6.011729	-108.440677	-4.225925	-1.0	True	1170.689441	0.163532	10975.160488	41.971349
297	35	FP all \| Scheurer	Impacted_Soil_Surface	Ocean_Surface_Water	6.018758	-108.696750	-4.448283	-1.0	True	1266.860208	0.753409	4509.321025	33.637734
298	4	FP all \| Scheurer	Impacted_Soil_Surface	Coast_Surface_Water	6.061087	-110.245151	-4.185854	-1.0	True	1183.945885	0.387262	4865.625691	162.674424
299	15	FP all \| Scheurer	Impacted_Soil_Surface	Coast_Surface_Water	6.093963	-111.455234	-4.489052	-1.0	True	1266.860208	0.753409	4509.321025	33.637734

	Case	Observed_dataset	Emission_Compartment	Target_Compartment	RMSE	R2	Slope	Spearman_r	Pass_Spearman	MPdensity_kg_m3	FI	t_half_deg_free	t_frag_gen_FreeSurfaceWater
220	6	Sea Surface all \| Isobe L2	Surface_Freshwater	Coast_Surface_Water	0.192740	0.540131	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
240	6	Sea Surface all \| Isobe L3	Surface_Freshwater	Coast_Surface_Water	0.196303	0.512105	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
121	137	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.210465	0.830557	-1.757989	-1.0	True	1052.671922	0.811533	3073.415520	14.554329
122	125	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.219214	0.816177	-1.605927	-1.0	True	1067.147583	0.814583	7649.971822	27.484251
260	26	Sea Surface all \| Isobe L4	Surface_Freshwater	Ocean_Surface_Water	0.219473	0.496743	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
200	26	Sea Surface all \| Isobe L1	Surface_Freshwater	Ocean_Surface_Water	0.219608	0.717863	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
123	135	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.220993	0.813182	-1.619336	-1.0	True	1266.860208	0.753409	4509.321025	33.637734
241	66	Sea Surface all \| Isobe L3	Coast_Surface_Water	Ocean_Surface_Water	0.225292	0.357365	-0.761238	-1.0	True	938.199911	0.224775	1827.090958	58.014555
124	121	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.226836	0.803171	-1.555892	-1.0	True	1027.531967	0.740818	4154.920691	47.156730
261	6	Sea Surface all \| Isobe L4	Surface_Freshwater	Coast_Surface_Water	0.228836	0.452885	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
125	131	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.240624	0.778517	-1.524742	-1.0	True	1193.888212	0.679468	7014.614289	66.504082
20	41	Sea Surface all \| Cai	Coast_Surface_Water	Coast_Surface_Water	0.241146	0.833343	-1.488332	-1.0	True	1027.531967	0.740818	4154.920691	47.156730
320	66	Sea Surface Fibers \| Song	Surface_Freshwater	Ocean_Surface_Water	0.242450	0.341708	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
140	123	Fibers in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.246943	-0.238814	-1.236710	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
221	66	Sea Surface all \| Isobe L2	Coast_Surface_Water	Ocean_Surface_Water	0.248879	0.233221	-0.761238	-1.0	True	938.199911	0.224775	1827.090958	58.014555
242	26	Sea Surface all \| Isobe L3	Surface_Freshwater	Ocean_Surface_Water	0.272301	0.061204	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
201	6	Sea Surface all \| Isobe L1	Surface_Freshwater	Coast_Surface_Water	0.273270	0.563134	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
321	46	Sea Surface Fibers \| Song	Surface_Freshwater	Coast_Surface_Water	0.274691	0.154986	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
80	46	Fragments in water \| Eo	Surface_Freshwater	Bulk_Freshwater	0.280510	0.837037	-1.644561	-1.0	True	938.199911	0.224775	1827.090958	58.014555

Evaluation of UTOPIA's performance in predicting MPs size dristributionsusing the monaco Monte Carlo framework¶

(https://monaco.readthedocs.io/en/latest/)¶

Import the necessary libraries¶

Load the UTOPIA's default configuration parameters¶

Define the parameter distributions for the Monte Carlo simulation.¶

Plot sensitivity indices for each output variable. Monaco uses the D-VARS approach to calculating global sensitivity indices.¶

Plot¶